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Abstract. We consider a delay differential equation (DDE) 
model for El-Nino Southern Oscillation (ENSO) variability. 
The model combines two key mechanisms that participate 
in ENSO dynamics: delayed negative feedback and seasonal 
forcing. We perform stability analyses of the model in the 
three-dimensional space of its physically relevant parame- 
ters. Our results illustrate the role of these three parameters: 
strength of seasonal forcing b, atmosphere-ocean coupling k, 
and propagation period r of oceanic waves across the Trop- 
ical Pacific. Two regimes of variability, stable and unstable, 
are separated by a sharp neutral curve in the (b, t) plane at 
constant k. The detailed structure of the neutral curve be- 
comes very irregular and possibly fractal, while individual 
trajectories within the unstable region become highly com- 
plex and possibly chaotic, as the atmosphere-ocean coupling 
k increases. In the unstable regime, spontaneous transitions 
occur in the mean "temperature" (i.e., thermocline depth), 
period, and extreme annual values, for purely periodic, sea- 
sonal forcing. The model reproduces the Devil's bleachers 
characterizing other ENSO models, such as nonlinear, cou- 
pled systems of partial differential equations; some of the 
features of this behavior have been documented in general 
circulation models, as well as in observations. We expect, 
therefore, similar behavior in much more detailed and realis- 



tic models, where it is harder to describe its causes as com- 
pletely. 

Keywords. Delay differential equations, El Nino, Extreme 
events, Fractal boundaries, Parametric instability. 

1 Introduction and motivation 

1 . 1 Key ingredients of ENSO theory 

The El-Nino/Southern-Oscillation (ENSO) phenomenon is 
the most prominent signal of seasonal-to-interannual climate 
variability. It was known for centuries to fishermen along 
the west coast of South America, who witnessed a seemingly 
sporadic and abrupt warming of the cold, nutrient-rich wa- 
ters that support the food chain in those regions; these warm- 
ings caused havoc to th eir fish harvests (Di az and Markgraf , 



1992; 



Philander 



19901) . The common occurrence of such 
warming shortly after Christmas inspired them to name it 
El Nino, after the "Christ child." Starting in the 1970s, El 
Nino's climatic effects were found to be far broader than just 
its manifestations off the s hores of Peru (IDiaz and Markgraf , 



1992; 



Glantz et al. 



1991). This realization led to a global 



awareness of ENSO's significance, and an impetus to at- 
tempt and i mprove prediction s of exceptionally strong El 
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Nino events (ILatif et al. 



199 A) 
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The following conceptual elements have been shown to 
play a determining role in the dynamics of the ENSO phe- 
nomenon. 



(i) The Bjerknes hypothesis: iBierknesI (11969b . who laid 
the foundation of modern ENSO research, suggested a pos- 
itive feedback as a mechanism for the growth of an internal 
instability that could produce large positive anomalies of sea 
surface temperatures (SSTs) in the eastern Tropical Pacific. 
We use here the climatological meaning of the term anomaly, 
i.e., the difference between an instantaneous (or short-term 
average) value and the normal (or long-term mean). Us- 
ing observations from the International Geophysical Year 
(1957-1958), Bjerknes realized that this mechanism must in- 
volve air-sea interaction in the tropics. The "chain reaction" 
starts with an initial warming of SSTs in the "cold tongue" 
that occupies the eastern part of the equatorial Pacific. This 
warming causes a weakening of the thermally direct Walker- 
cell circulation; this circulation involves air rising over the 
warmer SSTs near Indonesia and sinking over the colder 
SSTs near Peru. As the trade winds blowing from the east 
weaken and thus give way to westerly wind anomalies, the 
ensuing local changes in the ocean circulation encourage fur- 
ther SST increase. Thus the feedback loop is closed and fur- 
ther amplification of the instability is triggered. 

(ii) Delayed oceanic wave adjustments: Compensating 
for Bjerknes's positive feedback is a negative feedback in 
the system that allo ws a return to cold er conditions in the 



basin's eastern part (lSuarez and Schopf 



1988b . During the 



peak of the cold-tongue warming, called the warm or El Nino 
phase of ENSO, westerly wind anomalies prevail in the cen- 
tral part of the basin. As part of the ocean's adjustment to 
this atmospheric forcing, a Kelvin wave is set up in the tropi- 
cal wave guide and carries a warming signal eastward; this 
signal deepens the eastern-basin thermocline, which sepa- 
rates the warmer, well-mixed surface waters from the colder 
waters below, and thus contributes to the positive feedback 
described above. Concurrently, slower Rossby waves prop- 



agate westward, and are reflected at the basin's western 
boundary, giving rise therewith to an eastward-propagating 
Kelvin wave that has a cooling, thermocline-shoaling effect. 
Over time, the arrival of this signal erodes the warm event, 
ultimately causing a switch to a cold, La Nina phase. 

(Hi) Seasonal fo rc ing: A growing b ody o f work 



Ghil and Robertson, 



Jin et al. 



199 A 



1996; 



2000; 



Chan g et al. 



Tziperman et al. 



1994 



1994 



1995 



1995b points 



to resonances between the Pacific basin's intrinsic air-sea 
oscillator and the annual cycle as a possible cause for the 
tendency of warm events to peak in boreal winter, as well 
as for ENSO's intriguing mix of temporal regularities and 
irregularities. The mechanisms by which this interaction 
takes place are numerous and intricate and their relative 



importance is not y et fully understood (ITziperman et al 



1995 



Battisti, 



1988) 



1 .2 Formulation of DDE models 

Starting in the 1980s, the effects of delayed feedbacks and 
external forcing have been studied using the formalism 
of delay differential equatio ns (DDE) (see, inte r alia , 



Bhatta charrva and Ghill (Il982b : 



Ghil and Childress (1987) 



Hale ( 



1977b : 



Nussbaum 



for g eoscience applications, and 
(1998) for DDE theory). Several DDE systems have been 
suggested as toy models for ENSO variability. Battisti and 



Hirst (11989b have considered the linear autonomous DDE 

dT/dt = -aT(t-r)+T, a > 0, r > 0. (1) 

Here, T represents the sea-surface temperature (SST) aver- 
aged over the eastern equatorial Pacific. The first term on the 
right-hand side (rhs) of (fl]i mimics the negative feedback due 
to the oceanic waves, while the second te rm reflects Bjerk 



nes's positive feedback. As shown in dBattisti and Hirst , 



1989), Eq. ([TJ reproduces some of the main features of a 



fully nonlinear coupled at mosphe r e-ocean model of ENSO 



dynam ics in the tropics ((Battisti, 



1988 



Zebiak and Cane, 



1987). 
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Suarez and Schopf dl988l) and Battisti and Hirst (119891) 
also studied a nonlinear version of ((TJ, in which a cubic non- 
linearity is added to the rhs of the equation: 



dT/dt = -aT(t-r)+T- T 3 , 



(2) 



where < a < 1 and r > 0. This system has three steady 
states, obtained by finding the roots of the rhs: 

T =0, T ia = ±{l-a) 1 ' 2 . 

The so-called inner solution To is always unstable, while 
the outer solutions T± t 2 may be stable or unstable depend- 
ing on the parameters (a, r). If an outer steady state is un- 
stable, the system exhibits bounded oscillatory dynamics; in 



dSuarez and Sc hopf , 



1988) it was shown numerically that a 



typical period of such oscillatory solutions is about two times 
the delay r. 

The delay equation idea was very successful in explain- 
ing the periodic nature of ENSO events. Indeed, the delayed 
negative feedback does not let a solution fade away or blow 
up, as in the ordinary differential equation (ODE) case with 
r = 0, and thus creates an internal oscillator with period de- 
pending on the delay and particular form of the equation's 
rhs. DDE modeling has also emphasized the importance of 
nonlinear interactions in shaping the complex dynamics of 
the ENSO cycle. At the same time, many important details 
of ENSO variability still had to be explained. 

First, a delayed oscillator similar to (Q]) or (0 typically 
has periodic solutions with well-defined periods. However, 
the occurrence of ENSO events is irregular and can only be 
approximated very coarsely by a periodic function. Second, 
El-Nino events always peak during the Northern Hemisphere 
(boreal) winter, hence their name; such phase locking cannot 
be explained by a purely internal delayed oscillator. Third, 
the value of the period produced by the delay equations devi- 
ates significantly from actual ENSO interevent times of 2-7 
years. The delay r, which is the sum of the basin-transit 
times of the westward Rossby and eastward Kelvin waves, 



can be roughly estimated to lie in the range of 6-8 months. 
Accordingly, model (0 suggests a period of 1 .5-2 years, at 
most, for the repeating warm events; this is about half the 
dominant ENSO recurrence time. 



Tziperman et al. y,994) have demonstrated that these dis- 
crepancies can be removed — still within the DDE frame- 
work — by considering nonlinear interactions between the 
internal oscillator and the external periodic forcing by the 
seasonal cycle. These authors also introduced a more re- 
alistic type of nonlinear coupling between atmosphere and 
ocean to reflect the fact that the delayed negative feedback 
saturates as the absolute value of the key dependent variable 
T increases; note that in (Q]l the feedback is linearly propor- 
tional to the delayed state variable T{t — r). Munnich et al. 



(1991) studied an iterated-map model of ENSO and made a 
detailed comparison between cubic and sigmoid nonlineari- 
ties. A s a result, the sigmoid type of nonlinearity was cho- 



sen in ( Tzi perman et al 



1994), resulting in the periodically 



forced, nonlinear DDE 

dT/dt = -a tanh [nT(t - n)] 

+ f3 tanh [«T(t-T 2 )] +7 cos(2 7rt). 



(3) 



Model ([3]) was shown to have solutions that possess an in- 
teger period, are quasiperiodic, or exhibit chaotic behavior, 
depending on the parameter values. The increase of solu- 
tion complexity — from period one, to integer but higher 
period, and on — to quasiperiodicity and chaos — is caused 
by the increase of the atmosphere-ocea n coupling parame- 



ter k. The study (ITziperman et al. 



1994) also demonstrated 



that this forced DDE system exhibits period locking, when 
the external, "explicit" oscillator wins the competition with 
the internal, delayed one, causing the system to stick to an 
integer period; see also the more detailed analysis of phase 
locking in th e intermediate coupled model (ICM) of Jin et al. 



(1 1994 



1996) 



To summarize our motivation for the choice of a "toy 
model," work during the past 30 years has shown that ENSO 
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dynamics is governed, by and large, by the interplay of 
several nonlinear mechanisms that can be studied in sim- 
ple forced DDE models. Such models provide a convenient 
paradigm for explaining interannual ENSO variability and 
shed new light on its dynamical properties. So far, though, 
DDE model studies of ENSO have been limited to linear sta- 
bility analysis of steady-state solutions, which are not typical 
in forced systems, case studies of particular trajectories, or 
one-dimensional scenarios of transition to chaos, varying a 
single parameter while the others are kept fixed. A major ob- 
stacle for the complete bifurcation and sensitivity analysis of 
such DDE models lies in the complex nature of DDEs, whose 
numerical and analytical treatment is much harder than that 
of their ODE counterparts. 

In this work we take several steps toward a comprehensive 
analysis of DDE models relevant for ENSO phenomenol- 
ogy. In doing so, we also wish to illustrate the complexity 
of phase-parameter space structure for even such a simple 
model of climate dynamics. 

In Section [2] we formulate our DDE model, provide ba- 
sic theoretical results for this type of DDEs, present the nu- 
merical integration method used, and describe several so- 
lution types and their possible physical interpretation. In 
Section [3] we proceed to explore fully solution behavior 
over a broad range of the model's three most physically rel- 
evant parameters. We reproduce several dynamical solu- 
tion features and bifurcation scen arios previously reported 
in the literature for both simpler (Saunders and Ghil, 



and more detailed ( Jin et al. 



1994 



1998 



1995 



1994 



Ghil and Robertson, 



Diikstra , 



199 



200i 



6; 



2001) 



Tziperman et al. 



Neelin et al. 



1994 



2005) models, report new ones, and describe 
the corresponding three-dimensional (3-D) regime diagram. 
This 3-D regime diagram includes large regions of very 
smooth parameter dependence, as well as regions of very 
sensitive dependence on the parameters; the neutral surface 
separating simpler from more complex behavior exhibits rich 
and apparently fractal patterns. 



The discussion in Section |4] highlights the possibility of 
spontaneous, intrinsic transitions between the presence or 
absence of intraseasonal, higher-frequency fluctuations, as 
well as of interdecadal, lower-frequency variability. Such 
higher- and lower-frequency variability accompanies the sea- 
sonal and interannual oscillations that dominate our model 
solutions. This coexistence of variabilities on several time 
scales affects not only the mean and period of the solutions, 
but also the distribution of extreme warm and cold events. 

An illustration of the complexity "burst" caused by in- 
troducing a scalar delay in a simple ODE is given in Ap- 
pendix [A] Appendix [B] contains the proof of a key theoret- 
ical result presented in Sect. [2] while Appendix ICl provides 
details on numerical procedures for DDEs. 

2 Model and numerical integration method 

2. 1 Model formulation and parameters 

We consider a nonlinear DDE with additive, periodic forcing, 

dh(t) 



dt 



= — a tanh [nh(t — t)] + b cos(2n u>t); (4) 



here t > and the parameters a,K,r,b, and u) are all real and 
positive. Equation is a simplified one-delay vers ion of the 



two-delay model considered by Tziperman et al. d 19941) ; it 
mimics two mechanisms essential for ENSO variability: de- 
layed negative feedback via the function tanh(K z) and peri- 
odic external forcing. As we shall see, these two mechanisms 
suffice to generate very rich behavior that includes several 
important features of more detailed models and of observa- 



(Philander 


1990 


Bierknes, 


1969 


Neelin et al. . 


1994 


1998) 



is left for future work. 

The function h(t) in (|4|i represents the thermocline 
depth deviations from the annual mean in the Eastern 
Pacific; accordingly, it can also be roughly interpreted as 
the regional SST, since a deeper thermocline corresponds 
to less upwelling of cold waters, and hence higher SST, 
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and vice versa. The thermocline depth is affected by the 
wind-forced, eastward Kelvin and westward Rossby oceanic 
waves. The waves' delayed effects are modeled by the 
function tanh [nh(t — r)]; the delay r is due to the finite 
wave velocity and corresponds roughly to the combined 
basin-transit time of the Kelvin and Rossby waves. The 
particular form of the delayed nonlinearity plays a very 
imp ortant role in the behavior of a DDE model. Munnich et 



al. dl99ll) provide a physical justification for the monotone, 
sigmoid nonlinearity we adopt. The parameter n, which 
is the linear slope of tanh(/cz) at the origin, reflects the 
strength of the atmosphere-ocean coupling. The forcing 
term represents the seasonal cycle in the trade winds. 

The model © is fully determined by its five parameters: 
feedback delay r, atmosphere-ocean coupling strength k, 
feedback amplitude a, forcing frequency u, and forcing am- 
plitude b. By an appropriate rescaling of time t and depen- 
dent variable h, we let u> = 1 and a = 1. The other three 
parameters may vary, reflecting different physical conditions 
of ENSO evolution. We consider here the following ranges 
of these parameters: < t < 2 yr, < /t < oo, < 6 < oo. 

To completely specify the DDE model we need to pre- 
scribe some initi al "history," i.e. the behavior of h(t) on the 



interval [— r, 0) dHalel 11977b . In most of the numerical ex- 
periments below we assume h(t) = 1, — r < t < 0, i.e. 
we start with a warm year. Numerical experiments with al- 
ternative specifications of the initial history suggest that this 
choice does not affect our qualitative conclusions. 

2.2 Basic theoretical results 

To develop some intuition about the dynamics of Eq. ©, we 
consider two limiting cases. In the absence of the feedback, 
a = 0, the model becomes a simple ODE and hence has only 
a sinusoidal solution with period 1. One expects to observe 
the same behavior for b/a 3> 1. In the absence of forcing, 



b = 0, we obtain a well-studied DDE 

h(t) = -g [h(t - t)} , g(z) = tanh(« z). (5) 

The character of the solutions of this equation depends 
strongly on the delay r. For small delays, one expects to 
see behavior reminiscent of the corresponding ODE with 
zero delay; the general validity of such "small- delay expec- 



tations" is analyzed in detail in (iBodnan 120041) . For larger 
delays the nonlinear feedback might produce more com- 
plex dynamics. These heuristic intuitions happen to be true: 
Eq. © has an asymptotic solution that is identically zero for 
t < To, and admits periodic solutions with pe riod 4t for 



r > tq, where the critical delay is tq 



Nussbaum, 



1979; 



Chow and W alter 



7r/(2 k) dCao 



1996; 



19881) . In addition, it 



is known that the null solution is the only stable solution 
for t < tq. At t = tq the system undergoes a Hopf bi- 
furcation, and the trivial steady state transfers its stability 
to a periodic solution. Among other solutions, an impor- 
tant role is played by the so-called slow oscillating solutions, 
whose zeros are separated by a dista nce of at least the delay 



t. In particular, Chow and Walther (11988b showed that pe- 
riodic solutions with period 4 r and the symmetry condition 
—h(t) = h[t — 2t) are exponentially asymptotically stable; 
that is, any other solution will approach one of these solu- 
tions at an exponential rate. Moreover, for r > tq, these 



soluti ons may be the only stable ones (iChow and Walter 



19881) . 

These results can help explain the observed near- 
periodicity of ENSO variability. If one takes the delay 
t to equal approximately the transit time of the traveling 
ocean waves, namely 6 to 8 months, then the 4r internal 
period of the ENSO oscillator becomes 2-3 years. This 
remark, together with our further observations in Sect. 12.41 
provides a good justification for the observed quasi-biennial 



oscillation in Tro 


pical P 


(Diaz and Marks 


ral 




1992 


1995; 


Ghil et al. 


2002 


). 



Philander 



1990; 



Jiang et al. 
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Below we summarize basic theoretical results about 
Eg. dp i. Follow i ng the traditional framework (Hale, 



1977 



space X = 
h : [-r,0) 
as 



Nussbauml 119981) . we consider the Banach 
C([— t, 0),M) of continuous functions 
• R and define the norm for h 6 X 



|| h ||=SU P {|/l(t)|, t£ [-T.0)}, 

where | • | denotes the absolute value in M. For convenience, 
we reformulate the DDE initial-value problem (IVP) in its 
rescaled form: 

= - tanh [k hit -r)l +b cos(27ri), t > 0, (6) 

at 

h(t) = 0(f), f 6 [-r, 0), 0(f) G X. (7) 

Proposition 1 (Existence, uniqueness, continuous depen- 
dence) For any fixed positive triplet (r, K, 6), f/ze /VP ©-(121) 
has a unique solution h(t) on [0, oo). This solution depends 
continuously on the initial data (f>(t), delay r and the rhs of 
© considered as a continuous map f : [0, T) x X — > M, 
/or any finite T. 

Proof. See Appendix[B] 

From Proposition [T] it follows, in particular, that the sys- 
tem ©-(|7]i has a unique solution for all time, which depends 
continuously on the model parameters (r, k, b) for any fi- 
nite time. This result implies that any discontinuity in the 
solution profile as a function of the model parameters indi- 
cates the existence of an unstable solution that separates the 
attractor basins of two stable solutions. 

Our numerical experiments suggest, furthermore, that all 
stable solutions of ©-(|7]l are bounded and have an infinite 
number of zeros. 

2.3 Numerical integration 

The results in this study are based on numerical integration of 
the DDE ©-(17]). We emphasize that there are important dif- 
ferences between numerical integration of DDEs and ODEs. 



The first important difference is in the requisite initial data. 
The solution of a system of ODEs is determined by its value 
at the initial point t = to. When integrating a DDE, terms 
like h(t — r) may represent values of the solution at points 
prior to the initial point. Because of this, the given initial 
data must include not only h(to), but also a "history" of the 
values h(t) for all t prior to to m the interval that extends as 
far back at the largest delay. 

Another important issue in solving DDEs arises when a 
delayed value of the argument falls within the current inte- 
gration step. In order to avoid limiting the step size to be 
smaller than the smallest delay, or, alternatively, to avoid ex- 
trapolating the previous solution, an iterative procedure must 
be used to obtain successive approximations of the delayed 
solution that will yield a satisfactory local error estimate. The 
implementation of this iterative procedure affects profoundly 
the performance of a DDE solver. 

These and other specific features of DDE numeri- 
cal integration require development of special software 
and often the problem-specific modification of such 
software. We used here the Fortran 90/95 DDE solver 



dde_solver of Shampine and Thompson (2006), available 
at 



http://www.radford.edu/~thompson/ffddes/ This solver 



implements a (5,6) pair of conti nuously embedded, explicit 



Runge-Kutta-Sarafyan methods (ICorwin et a 



1997). Tech- 



nical details of dde_solver, as well as a brief overview of 
other available DDE solvers are given in Appendix [C] 

The numerical simulations in this paper require very long 
integration intervals, leading to prohibitive storage require- 
ments. This difficulty led us to incorporate several new 
options in dde_solver; they are also described in Ap- 
pendixlcl 

2.4 Examples of model dynamics 

In this subsection we illustrate typical solutions of the prob- 
lem ©-(17]) and comment on physically relevant aspects of 
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Fig. 1. Examples of DDE model solutions. Model parameters are 
k = 100 and 6=1, while k increases from curve (a) to curve (f) as 
follows: (a) r = 0.01, (b) r = 0.025, (c) r = 0.15, (d) r = 0.45, 
(e) t = 0.995, and (f) r = 1. 

these solutions. All experiments shown here use the constant 
initial data 0=1. Figure Q] shows six trajectories obtained 
by fixing b = 1,k — 100 and varying the delay r over two 
orders of magnitude, from r = 10~ 2 to t = 1, with r in- 
creasing from bottom to top in the figure. The sequence of 
changes in solution type illustrated in Fig.[T]is typical for any 
choice of (b, k) as r increases. 

For a small delay, r < 7r/(2 k), we have a periodic so- 
lution with period 1 (curve a); here the internal oscillator 
is completely dominated by the seasonal forcing. When the 
delay increases, the effect of the internal oscillator becomes 
visible: small wiggles, in the form of amplitude-modulated 
oscillations with a period of 4r, emerge as the trajectory 
crosses the zero line. However, these wiggles do not affect 
the overall period, which is still unity. The wiggle ampli- 
tude grows with r (curve b) and eventually wins over the 
seasonal oscillations, resulting in period doubling (curve c). 
Further increase of r results in the model passing through 
a sequence of bifurcations that produce solution behavior of 
considerable interest for understanding ENSO variability. 



Some of these types of solution behavior are illustrated 
further in Fig. [2] Panel (a) (k = 5, r = 0.65) shows the 
occurrence of "low-/i," or cold, events every fourth seasonal 
cycle. Note that negative values of h correspond to the boreal 
(Northern Hemisphere) winter, that is to the upwelling sea- 
son — December-January-February — in the eastern Trop- 
ical Pacific; in the present, highly idealized model, we can 
associate the extreme negative values with large-amplitude 
cold events, or La Ninas. This solution pattern loses its regu- 
larity when the atmosphere-ocean coupling increases: Panel 
b (k = 100, t = 0.58) shows irregular occurrence of large 
cold events with the interevent time varying from 3 to 7 cy- 
cles. 

In panel c (k = 50, r = 0.42) we observe alter- 
nately and irregularly occurring warm El-Nino and cold 
La Nina events: the "high-/i" events occur with a pe- 
riod of about 4 years and random magnitude. Panel d 
(k = 500, t — 0.005) shows another interesting type 
of behavior: bursts of intraseasonal oscillations of ran- 
dom amplitude superimposed on regular, period-one dynam- 



tions ( 


Madden and Julian, 


1971, 


1972 


1994 


) or westerly- 


wind bursts 




Gebbie et al. 


2007; 


Harrison and Giese 


1988; 


Verbickas, 


1998 




Delcroix et al. 


1993 


). The solution in 



panel e (k = 50, r = 0.508) demonstrates sustained inter- 
decadal variability in the absence of any external source of 
such variability. The solution pattern illustrates spontaneous 
changes in the long-term annual mean, as well as in the dis- 
tribution of positive and negative extremes, with respect to 
both time and amplitude. 

3 Critical transitions 

3.1 Numerical characterization of solution behavior 

In this section we focus on the onset of instabilities in 
the model (|6]l-(|7]). Taking a "metric" approach to the 
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Fig. 2. Noteworthy solution patterns of relevance to ENSO dynam- 
ics; seasonal forcing amplitude b — 1. a) Regularly occurring cold 
(low-/i) events, or La Ninas (k — 5, t — 0.65); b) irregular cold 
events (ft = 100, r = 0.58); c) irregular alternations of warm (El 
Nino, high-fr) and cold events (ft = 50, r = 0.42); d) intraseasonal 
activity reminiscent of Madden- Julian oscillations or westerly-wind 
bursts (ft = 500, r = 0.005); and (e) interdecadal variability in the 
annual mean and in the relative amplitude of warm and cold events 
(ft = 50, r = 0.508). 

problem, we study the change in several statistics of a 
trajectory as the model parameters change. This approach 
is complementary to the "topological" one, which forms 
the basis for the stability a nalysi s of dynamical systems 



( Andronov and Pontrvagin , 



1937 



Katok and H asselblatt. 



1995). In the latter, one studies the topological structure 
of the system's attractor(s), i.e. a combination of points, 
circles, tori, or more complicated objects. The motivation 
for this approach comes from noting that topologically 
equivalent solutions can be mapped onto each other using an 
appropriate diffeomorphism, i.e. a one-to-one, continuously 
differentiable map. Hence, considering topologically equiv- 
alent classes rather than all individual solutions is enough for 



studying the system's qualitative behavior. From a practical 
point of view, though, metric properties of a solution might 
be as important as its topological description or more: think 
of living in a region with constant air temperature of -10°C 
vs. 20°C. Furthermore, metric properties are also much 
easier to study, in simple models, as well as in full, 3-D 
general circulation models (GCMs) and in observational 
data sets. 

Technically, we proceed in the following way. For each 
fixed triplet of parameters (b, k, t) we find a numerical ap- 
proximation hi = h(ti) of the model solution h(t) on a 
grid Q = {tj}i=i,...,jv, U € [0, T max ] using the Fortran 
90 DD E solver dde_solver of Shampine and Thompson 
(2006) (see also Appendix P. We only consider the latter 
part of each solution, ti > T m - m > 0, in order to avoid any 
transient effects; to simplify notations, we assume from now 
on that t\ = T m i n . Typical parameters for our numerical 
experiments are T max = 10 4 , T max — T m ; n = 10 3 , time 
step 5 = ti — ti-x = 10 -3 , and a numerical precision of 
e = 10~ 4 . We have also verified some of the results with a 
precision up to 10 -12 , a time step of 10~ 4 , and over time in- 
tervals up to T max — T m i n = 10 4 in order to ascertain that the 
reported phenomena are not caused or affected by numerical 
errors. 

We report results for the following trajectory statistics: 
maximum value M = maxj hi, mean value E = J^i hi/N, 
where N is the maximal integer less than (T max — T m i n )/S; 
and mean of positive values E + = J2i W/^X)/ Si 1 i H >a' 
where 1^ is the characteristic function of the set A, identi- 
cally equal to one for all points in A and zero outside of A. 
Furthermore, we have computed, but do not show here, the 
trajectory variance, the mean of negative values, and upper 
90% and 95% quantiles; the results are very similar for all 
the statistics we have examined. 

We also computed an approximation to the period of a so- 
lution. Specifically, for any positive integer A we define the 
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A-discrepancy 

E^Ia+i (hi - hi-A) 

Ra = (8) 

(N - A) Var(h) 

1 N 2 1 N 

Var{h) = -Y: (hi-h) , h=- y £h i . 

i=i i=i 

If hi = h(U), for a periodic solution h with period P = S A, 
we have Ra = 0. In numerical experiments, we can only 
guarantee that the A-discrepancy of a periodic solution is 
small enough: Ra < r e — 4e 2 /Var(h), where e is the ab- 
solute numerical precision. We call near-period the minimal 
number V = S A such that 

Ra < r r] and Ra < Ra±i 

for some prescribed < rj -C 1. The first condition ensures 
that the A-discrepancy Ra is small enough, while the second 
one guarantees that Ra is a local minimum as a function 
of A. The following proposition follows readily from the 
definition of V. 

Proposition 2 (Convergence theorem). The near-period V 
converges to the actual period T of a continuous periodic 
function h(t) when the numerical step 5 and nominal accu- 
racy e decrease: lim^o. e^o V = T. 

The continuity requirement in the hypothesis is too re- 
strictive for this statement, but it suffices for our purpose, 
since solutions of model (O, for a given triplet (b, t, k), are 
smooth. The rate of convergence in Proposition [2] depends 
strongly on the period structure of h(t). The rate is high for 
functions with "simple" periods, e.g. with a single local max- 
imum within the period T, and may be arbitrarily low in the 
general case; e.g. for h = hi + 7/12, where hi has period 
T/2, /12 period T, and 7 is small enough, the convergence is 
quite slow. 

To summarize, the near-period V approximates the actual 
period T for periodic functions. The functional V can also 
be defined for certain functions that are not periodic. For 
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Fig. 3. Near-period V as a function of delay r for fixed b = 0.03, 
heavy solid segments; dashed straight lines correspond to V = 
Ark, k — 1,2, .... The near-period is always a multiple of 4r 
close to an integer. 

instance, this may be the case for a quasiperiodic function 
h(t) — pi(t)+p2(t), where each pi (t) is periodic and the two 
periods, T\ and T2, are rationally independent. It is also the 
case for a near-periodic function h(t) = Po(t) +Pi(t), where 
Po (t) is periodic, and p\ (t) has sufficiently small amplitude. 
As we shall see, the period approximation V is quite helpful 
in understanding the structure of our model's solution set. 

3.2 Small forcing amplitude and frequency locking 

We mentioned in Sect. 12.21 that, without external forcing 
(b = 0), the nontrivial stable solutions of the model ©-(|7]i 
are periodic with period 4 r. When the external forcing is 
small, b <C 1, the dynamical system tries to retain this prop- 
erty. Figure [3] shows the near-period V as a function of the 
delay r for fixed b = 0.03 and k — 100. Here r v = Ri, 
that is we compare the A-discrepancy Ra with the one-step 
discrepancy R\ ; the latter measures the degree of continuity 
of our discrete-time approximation to h(t). Straight dashed 
lines in the figure correspond to V — 4 r k for positive inte- 
ger k. One can see that the solution's near-period is always a 
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3.3 Onset of the instabilities 




0.005 0.01 0.O15 0.02 0.025 0.03 
Forcing amplitude, b 



Fig. 4. Devil's bleachers: period index k = V /{At) as a function 
of forcing amplitude b and delay r. Notice the presence of very 
long periods, of 100 years and more, in the color bar and figure. 



multiple of 4 r and always close to an integer. This state of 
affairs is a natural compromise between the internal period 
At and the driving period 1, a compromise rendered possi- 
ble by the internal oscillator's nonlinearity. 

More generally, Fig. [4] shows the map of the period 
index k = P/(At) in the b-r plane for < b < 0.03, 
< r < 0.22. Here one immediately recognizes 
the so-called Devil's bleachers scenario of transition 
to chaos, documented in other ENSO models, includ- 



ing b o th ICM s ( Jin et al 



199 A 



1995 



1994 



Ghil and Robertson, 



Ghil and Robertson, 



1996; 



Tziperman et al. 



2000) and GCMs 



2000), as well as in certain obser- 



vations dGhil and Robertson , 



2000; 



Yanai and Li , 



1994). 



The periodically forced model (|6]l exhibits the web of 
resonances that characterizes coupled oscillators, although 
the "external oscillator" quickly wins over the internal one 
as the amplitude b of the forcing increases. For b = 0.1 (not 
shown) the near-period plot looks similar to the one shown 
in Fig. [3] but V takes only integer values. When the forcing 
amplitude increases further, the near-period V, as well as the 
actual period T, is locked solely to integer values. 



Tziperman et al. (1994) reported that the onset of chaotic 
behavior in their two-delay, periodically forced DDE model 
is associated with the increase of the atmosphere-ocean cou- 



pling k; Munnich et al. (119911) made a similar observation 



for an iterated-map model of ENSO. We explore this tran- 
sition to chaos in our model over its entire, 3-D parameter 
space. 

First, we compute in Fig.|5]the trajectory maximum M as 
a function of the parameters b and r for increasing values 
of K. For small values of coupling (top panel) we have a 
smooth map, monotonously increasing in b and periodic with 
period 1 in r. As the coupling increases, the map loses its 
monotonicity in b and periodicity in t for large values of t, 
but it is still smooth. For k w 2 (middle panel), a neutral 
curve f(b, r) = emerges that separates a smooth region (to 
the right of the curve), where we still observe monotonicity 
in b and periodicity in r, from a region with rough behavior 
of M. The gradient of M(b, r) is quite sharp across this 
neutral curve. 

Further increase of the coupling results in a qualitative 
change in the maximum map. The neutral curve, which be- 
comes sharp and rough, separates two regions with very dif- 
ferent behavior of M(b, r) (bottom panel). To the right of 
the curve, the map M(b, r) is still smooth, periodic in r and 
monotonic in b. To the left, one sees discontinuities that pro- 
duce rough and complicated patterns. The mean position of 
the neutral curve f(b,r) = quickly converges to a fixed 
profile, although its detailed shape at smaller scales contin- 
ues to change with increasing k. The limiting profile is close 
to the one observed for n = 11 (bottom panel). 

3.4 Unstable behavior 

In this subsection we illustrate the model's parametric in- 
stabilities using the four trajectory statistics introduced in 
Sect. 13.11 maximum M, mean E, mean of positive values 
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a) k=0 5 




n ? 4 b a m 



Forcing amplitude, b 
b) w=2 ,0 




Forcing amplitude, b 



c)k=11.0 




2 4 6 3 10 

Forcing amplitude, b 



Fig. 5. Maximum map M = M(b, r). Top: n = 0.5, middle: 
k = 2, and bottom: « = 11. Notice the onset of instabilities and 
emergence of a neutral curve f(b, r) = that separates the smooth 
from the unstable regions. 
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E+, and near-period V . Figure [6] shows a plot of these statis- 
tics in a rectangle of the plane (b, r) for fixed k — 10. The 
neutral curve of Fig. [5J; crosses this rectangle from its bot- 
tom left corner to the central point on its right edge; thus 
the bottom right region of each panel corresponds to smooth 
behavior of each statistic map, while the top left region cor- 
responds to rough behavior. This figure illustrates the fol- 
lowing points: 

The maximum map M(b, r) (Fig. [(J) shows, among other 
instabilities, a pronounced jump along a mainly horizontal 
curve in the (b, r) plane. In the vicinity of this curve, an ar- 
bitrarily small increase in r causes a 200-300% jump (from 
0.25 to 0.5-0.8) in M, A pair of trajectories on either side of 
this transition is shown below in Fig. [9] 

The mean map E(b, r) (Fig. |6j?). To emphasize the para- 
metric instabilities, we show here log 10 \E\. Deviations of 
E from reflect the trajectory's asymmetry; hence the larger 
values of this map indicate asymmetric solutions. In this ex- 
periment, we use a numerical precision of 10~ 4 , so that val- 
ues of log 10 \E\ < —4 effectively correspond to symmetric 
trajectories, E = 0. One can see that the symmetry, charac- 
teristic for trajectories from the smooth region (bottom right 
part), breaks across the neutral curve. In the unstable re- 
gion, the magnitude of the asymmetry is very intermittent; 
it ranges over three orders of magnitude, taking its maximal 
value in the region that corresponds to the jump in the trajec- 
tory maximum, cf. panel (a). 

The mean of positive values E + (b, r) (Fig. [6};), by com- 
parison with the maximum map in panel (a), shows that cer- 
tain internal instabilities may affect trajectory shape without 
affecting the behavior of extremes. For instance, the maxi- 
mum map M (b, r) is smooth within the neighborhood of the 
point (b = 1.5, r = 0.49), although the map E + (b, r) ex- 
hibits a discontinuity across this neighborhood. In fact, one 
arrives at the same conclusion by comparing the maximum 
map to the mean map, cf. panels (a) and (b), respectively. 
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b) Mean, E 
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Forcing amplitude, b 
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c) Mean of positive values, £ + 
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d) Near-period, P 
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Fig. 6. Trajectory statistics plots for k = 10, as a function of forcing amplitude b and delay r: a) maximum map, M(b, r); b) mean map, 
log \E(b, t)|; c) mean of positive values, E+(b, r); and d) near-period map, V(b, r). 



The near-period map V(b, t) (Fig.|6}J). The near-period 
is varying over the interval [0, 27] in this map. As we have 
noticed, not all of these values correspond to trajectories that 
are actually periodic, rather than just nearly so (see Sect. 13. Il l 
. The large constant regions, though, do reflect the actual 
periods; as a rule, they correspond to small values of V. Ex- 
amples include: V = 1, within the smooth part of the map 
(bottom right); V = 2 within the middle horizontal tongue; 
V = 3 within the top right part; and V = 5 in a small tongue 
that touches the left margin of the plot at (b = 1, r = 0.44). 

Figure [7] shows the same plots of solution statistics over a 
rectangle of the plane (k, t), for a fixed value b = 1; it il- 



lustrates the onset and development of instabilities as a func- 
tion of the coupling parameter k. Comments similar to those 
above, which concern the behavior of individual maps and 
the connections between them, apply to this set of maps as 
well. 

Figure[8]shows the maximum and mean maps for n = 100. 
The increase in the atmosphere-ocean coupling, with its as- 
sociated nonlinearity, amplifies the model's instabilities and 
leads to more complex dynamics that is quite chaotic. The 
rigorous verification of the chaotic properties is left, how- 
ever, for future work. 

Figure [9] shows three examples of change in solution be- 
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a) Maximum, M 
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Atmosphere-ocean coupling, k 



c) Mean of positive values, £ + 
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Atmosphere-ocean coupling, k 
Fig. 7. Solution statistics plots for b = 1, as a function of couplin; 



b) Mean, E 




5 10 15 20 25 30 35 40 45 



Atmosphere-ocean coupling, k 



d) Near-period, P 




5 10 15 20 25 30 35 40 45 



Atmosphere-ocean coupling, re 

parameter k and delay r. Same panels as in Fig. [6] 



havior across the neutral curve that separates smooth from 
rougher behavior. In all three cases shown, one trajectory 
(dashed line) has period one and lies within the smooth part 
of the parameter space, in the immediate vicinity of the neu- 
tral curve, while the other trajectory (solid line) corresponds 
to a point in parameter space that is quite close to the first 
one but on the other side of the neutral curve. 

Panel (a) illustrates transition to quasiperiodic behavior 
with a "carrier wave" of period near 8, alternating smoothly 3 
warm and 3 cold events, separated by one "normal" year. In 
panel (b) we see single large El Ninos alternating with sin- 
gle large La Ninas, separated by 3 normal years. Panel (c) 
exhibits interdecadal variability, like in Fig. 2e, except that 



the transitions between warm and cold "decades" is sharper 
here, and the spells of El Ninos and La Ninas even longer (a 
dozen years here vs. roughly ten there). 

4 Discussion 



We have considered a toy model for ENSO variability that 
is governed by a delay differential equation (DDE) with a 
single, fixed delay and periodic forcing. Thus, we follow 
a line of research pioneered by Suarez and Schopf (1988), 
Battisti and Hirst Jl989h . and Tziperman et al. i 1994 . who 
have shown that DDE models can effectively capture com- 
plex phenomena found in much more detailed ENSO mod- 
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Fig. 8. Solution statistics plots for k — 100: a) maximum map, M(b, r); and b) mean map, log \E(b, r)|. 



els, as well as in observational data sets. DDE models are 
very simple and, at the same time, exhibit rich and complex 
behavior. Stability and bifurcation analysis for such models 
can be carried out analytically only to some e xtent, but nu 



merical methods a r e being actively developed (iBaken. l2000t observations (IJin et al 



Baker et al. 



1995 



for ENSO dynamics, we emphasize the following observa- 
tions. A simple DDE model ©-(|7]) with a single delay re- 
produces the Devil's staircase scenario documented in other 
ENSO model s, includ i ng ICMs and GCMs, as wel l as in 



Engelborghs et al. 



2001), and we have 



not yet taken full advantage here of either approach. 

To initiate stability and bifurcation analysis of ENSO- 
related DDE models, we started here with a descriptive nu- 
merical exploration of Eqs. ©-(IT]) over a wide range of phys- 
ically relevant parameter values. We studied parameter de- 
pendence of various trajectory statistics, and report the exis- 
tence of a large domain in parameter space where the statis- 
tics maps are strikingly discontinuous. The local continuous- 
dependence theorem (Proposition [T]i suggests, at least, that 
the reported discontinuities in global solutions point to the 
existence of unstable solutions of Eqs. (|6]l-(|7]l; the complex 
discontinuity patterns (see Figs. [6] and [7]) lead us to suspect 
the presence of a rich family of unstable solutions that under- 
lie a complicated attractor. It would be interesting to study in 
greater detail the unstable solutions and understand their role 
in shaping the system dynamics. 

Summarizing the model results in terms of their relevance 



1995 



1994 



Ghil and Robertson, 



1996 



Tziperman et al. 



199 A 



2000). The model illustrates, in 



simplest possible terms, the role of the distinct parameters: 
strength of seasonal forcing b vs. ocean-atmosphere cou- 
pling k and transit time r of oceanic waves across the Trop- 
ical Pacific. We find spontaneous transitions in mean ther- 
mocline depth, and hence in sea surface temperature (SST), 
as well as in extreme annual values that occur for purely 
periodic, seasonal forcing. The model generates intrasea- 
sonal oscillations of various periods and amplitudes, as well 
as interdecadal variability. The f ormer result might suggest 



that Madde n-Julian oscillations (Madden and Julian, 



1972, 



1994) and w e sterly wind burst s 



Harr: 



son and Giese 



1988 



Verbickas, 



Gebbie et al 



1998 



1971 



2007 



Delcroix et al. 



19931) in the Tropical Pacific are affected by ENSO's interan- 
nual modes at least as much as they affect them in turn. The 
latter result likewise suggests that interdeca dal variability in 



the extratropical, fherm ohaline circulation (iDijkstra . 



2005 



Dijkstra and Ghil , 



20051) might also interfere constructively 
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Fig. 9. Examples of instabilities across the neutral curve that sepa- 
rates smooth from rough behavior; at b = 1; see Fig. [5] The dashed 
blue curve shows a period- 1 trajectory on the smooth side of the 
neutral curve; the solid red curve shows a trajectory immediately 
across the neutral curve, in the rough-solution domain. 

with ENSO's intrinsic variability on this time scale. 

A sharp neutral curve in the (b — r) plane separates smooth 
parameter dependence in the solution s' map of "metrics" 



(Taylor 



2001 



Fuglestvedt et al. 



20031) from "rough" behav- 



ior. We expect such separation between regions of smooth 
and rough dependence of solution metrics on parameters in 
much more detailed and realistic models, where it is harder 
to describe its causes as completely. 

Finally, it appears that, even in as simple a model as 



our DDE, the mean, extrema and periodicity of solutions 
can change (a) spontaneously, without any change in the 
external forcing; and (b) one of these characteristics can 
change considerably, while others change but very little. 
Furthermore, certain parts of parameter space involve only 
small and smooth changes, while others involve large and 
sudden ones. It is quite conceivable that such behavior 



might arise in intermedia te climate models (Jin et a 



Neelin et al. 



1994 



Stainforth et al. 



2005). 



1998b and GCMs dMurphv et al. 



1996; 



2004 



Appendix A 



A simple example of DDE complexity 



In his classical book (Hale, 



1977b on functional differential 



equations, Jack Hale remarks that systematic study of dif- 
ferential equations with dependence on the past started with 
the work of Volterra on predator-prey models and viscoelas- 
ticity at the beginning of the 20th century. DDEs have thus 
been actively studied and applied for almost a century. Still, 
they are a relatively new modeling tool when compared to 
ODEs, and their theory and numerical analysis are much less 
developed than for ODEs. To develop the reader's intuition 
for DDEs, we discuss in this appendix a simple autonomous 
ODE and the corresponding DDE obtained by introducing a 
fixed time delay; our goal is to illustrate how this apparently 
innocuous modification complicates the solution set of the 
equation and renders its analytical and numerical study more 
elaborate. 

We start with the linear, scalar ODE 



x(t) = a x(t). 



(Al) 



Assuming a solution of the form x — ce At , we substitute it 
in (IAU to find its characteristic equation 

A = a. 
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This equation allows us to find all possible functions of the 
given form that satisfy ( lAU . Clearly, fixing an initial con- 
dition x(Q) = xq leaves only one solution with c = xq. 
This example illustrates an important general property of au- 
tonomous ODEs: their characteristic equations are polyno- 
mials in A and thus have a finite set of (complex-valued) so- 
lutions that can be easily found. As a result, the finite set of 
solutions to (lAlb can also be described easily. 
Let us introduce now a delay r into Eq. (IA1I) : 



x(t) = a x(t — t). 



(A2) 



This modification implies that it takes some finite time r for 
changes in the model state x(t) to affect its rate of change 
x(t). Such an assumption makes sense in many applica- 
tions, with nume rous sp ecific examples given in the Intro- 



duction to Hale dl977l) . and in Kolmanovskii and Nosov 
(1986). Proceeding as before, we assume a solution of the 
form x = c e A * and obtain the characteristic equation 



(A3) 



The fact that not all exponential terms cancel out in ( IA3b 
changes dramatically the solution set of this characteristic 
equation. It can be shown that dA3b has an infinite number of 
complex solutions; hence there exists an infinite numbe r of 



functions that satisfy Eq. (1A2I) (Hale, 



general solution to (IA2b is given by ( Falb 



1977 



Falbo 



(l ll995h : 



1995 ). A 



x(t) = Co e~ % l r + Ci e r " 1 + C 2 e ri * + C 3 e r 1 

OO 

+ eP " 1 [^1* cos (* *) + °2k sm(q k t)] 



(A4) 



fc=i 



Here pk ± iqk are complex solutions of ( lA31 l. C\k and C^k 
are arbitrary constants, and Co, Ci, 62,63 depend on the 
values of a and r: for a < — l/(re), C, = 0, i = 0, 1, 2, 3; 
for a — — l/(re), Cq is arbitrary, and C,- = 0, i = 1, 2, 3; 
for -l/(re) < a < 0, C = C 3 = 0, Ci,C 2 are arbitrary, 
and ro, ri are the real roots of ( 1A3I ): finally for a > 0, Ci — 
for i = 0, 1, 2, C3 is arbitrary, and r is the only real root of 



dA3b . Accordingly, this simple autonomous linear DDE may 
have increasing, decreasing, and oscillating solutions. 

The "burst" of complexity in the solution set of (lA2b com- 
pared to that of (IA1I) is pervasive in the DDE realm. Solving 
characteristic equations arising from DDEs is typically quite 
involved and requires usually some general results on equa- 
tions that mix polynomial and exponential terms (see, e.g., 



Appendix A in 



Hald (11977J)). Stability and bifurcation stud- 



ies for DDEs face a similar "burst" of complexity, compared 
to ODEs. 

Numerical exploration of DDEs is also considerably more 
complicated than for ODEs. The question of discontinuities 
in the solutions' derivatives is crucial: Because general nu- 
merical methods for both ODEs and DDEs are intended for 
problems with solutions that have several continuous deriva- 
tives, discontinuities in low-order derivatives (up to the order 
of the integration method) require special attention. Such 
discontinuities are not all that rare for ODEs, but they are al- 
most always present for DDEs, since the first derivative of the 
history function at the initial point is almost always different 
from the first derivative of the solution there. Moreover, dis- 
continuities are a much more serious matter for DDEs than 
they are for ODEs because they propagate: the smoothness of 
the derivative x at the current time t depends on the smooth- 
ness of the solution x at past times in the interval t — r. 

This difficulty may be illustrated using our very simple 
example 

x(t)=x(t-l), (A5) 

where we set, without loss of generality, r = a = 1. For 
this equation x^ k+1 ^(t) — x^ k \t — 1), where x^ denotes 
the k-t\\ derivative. In general, if there is a discontinuity of 
order k at a time t*, meaning that x^ has a jump at t = t*, 
then as t crosses the value t* + 1, there is a discontinuity in 
^Ct'+i) because of the term x(t — 1). With multiple delays r,, 
a discontinuity at time t* propagates to the times 



t*+n, t* 



T2, 



t*+T k , 
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and each of these discontinuities in turn propagates. If there 
is a discontinuity of order k at time t*, the discontinuity at 
each of the times t* + Tj is at least of order k + 1. Because 
the effect of a delay appears in a higher-order derivative, the 
solution does become smoother as the integration proceeds. 
This "generalized smoothing" proves to be quite important 
to the numerical solution of DDEs. 

To illustrate these statements, suppose we wish to solve 
Eq. ( IA5t with history x(t) = 1 for t < 0. On the interval 
< t < 1 the solution is x(t — 1) = 1 because t — 1 < 0. 
Thus, the DDE on this interval reduces to the ODE x(t) = 1 
with initial value x(0) = 1. We solve this problem to obtain 
x(t) = t+ 1 for < t < 1. Notice that this solution exhibits 
a typical discontinuity in its first derivative at t = because 
it is zero to the left of the origin and unity to the right. Now 
that we know the solution for t < 1, we can reduce the DDE 
on the interval 1 < t < 2 to an ODE x = (t - 1) + 1 = t 
with initial value a;(l) = 2 and solve this problem to find 
that x(t) = Q.5t 2 + 1.5 on this interval. The first derivative is 
continuous at t = 1, but there is a discontinuity in the second 
derivative. In general, the solution of the DDE on the interval 
[k, k + 1] is a polynomial of degree k + 1 and the solution 
has a discontinuity of order k + 1 at time t = k. 

In order for a solver to account for these specific features 
of DDEs and to solve them efficiently, accurately and reli- 
ably, there must be a great deal of care taken "under the 
hood" of the solver. Discontinuities need be tracked only 
up to the order of the integration method, since higher-order 
discontinuities do not affect the performance of the solver, 
i.e. its procedures for error estimation and step size control. 
Our solver of choice, dde_solver, tracks discontinuities 
explicitly and includes them as integration grid points, in or- 
der to avoid interpolating across them. For problems with 
constant delays, it is possible to build the necessary "discon- 
tinuity tree" in advance and to step exactly to each point in 
the tree. 

For problems with state-dependent delays, the discon- 



tinuity times are not known in advance. In this case, 
dde_solver tracks discontinuities using root finding, in 
conjunction with the primary integration method's under- 
lying polynomial interpolants, to locate the discontinuities 
and restart the integration at each such point. Some other 
available solvers handle discontinuity propagation differ- 
ently; but the best solvers do take special precautions of 
one kind or another, since ignoring discontinuities can sig- 
nificantly affect the reliability of a DDE solver. We refer 
in this context also to several distinct approa ches to the nu- 

equations i 



1984: 


Ghil and Mullhauot 


Zaliapin et al. 


2003 





1985 



Dee and Ghil! 



S au nders and G hil , 



2001; 



Appendix B 

Proof of Proposition Q] 

Consider the IVP ©-([7]), with the rhs of the DDE denoted 
by 

T{t,h(-)) = -t&nh[Kh(-)] +b cos(2 7ri). 

Existence of the solution to this problem on [0, T] for some 
T > readily follows from the continuity of F(t,h ) and 



Hale 



1977), The- 



the general existence theorem for DDEs 
orem2.1,p. 41). Moreover. Nussbaum dNussbauml (11998). 
Theorem 1, p. 3) remarks that if there exist constants A and 
B such that || T(t,h) ||< A || h || +B for all (t, h) elxl 
then one can choose T = oo. Since | tanh(Kz)| < k\z\, the 
choice of A = k and B — b ensures the solution's existence 
on [0, oo). 

Uniqueness could be derived from the Lipschitz pro perty 
of F{ t, h) in h and the general uniqueness theorem (Hale 
( 1977), Theorem 2.3, p. 42). However, for our system the 
uniqueness can be established in a simpler way. Indeed, as- 
sume that x(t) and y(t) are solutions of ©-(O on [0, T], 
with rhs JF(i, h(t — t)) and the initial condition </>(<). Then, 
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for < t < t, 
x (t) ~ y( f ) = I [F{s, x(s - t) - F(s, y(s - t)] ds 



(5,6) pair of continuously embed ded, explicit Runge-Kutta 



(tanh[Ky(s — r)] — tanh[/ix(s — t)]J ds 

— T 

( tannic — tanh[«;x(it)]) du = 0. 



(Bl) 



Thus, the solutions x and y are identical up to t = r > 0. 
The uniqueness is proven by successively advancing in time 
by intervals of length r. 

Continuous dependence on initial conditions and the rhs 
of (O for any finite T follows from the existence and un ique 



ness a nd the general continuous dependence theorem (Hale 



(1977), Theorem 2.2, p. 41). To show the continuous de- 
pendence on the delay r, also for finite T, we consider the 
sequence r k — > t, k — » oo, and for any r k introduce the time 
scale change: s k = tr/rk. Then, 

h(t - Tk) = h (-^( s k - t)j = hk{sk - t), 



and one finds 

_d_ 

dsk 



hk(sk) = ~—A K [h k (s k -r)] + b cos (2?r^-^ 



Clearly, the rhs of this system converges to T as k — > cxo. 
This shows that a small change in the delay r can be consid- 
ered as a small change of the rhs of the Eq. (O with the same 
delay. Hence, continuous dependence on the rhs implies con- 
tinuous dependence on the delay. 

Appendix C 

DDE solvers 

CI Our solver of choice: dde_solver 



The D DE solver dde_solver (IShampine and Thompson , 
2006) was used to perform the numerical experiments re- 
ported in this paper; dde_solver is a Fortra n 90/95 exten- 



sion o f its Fortran 77 predecessor dklag6 (ICorwin et al 



Sarafyan methods. We refer to (ICorwin et al. 



1997) for the 



coefficients and precise details of the methods used, and to 



dShampinei 1 1 9941) for a discussion of continuously embedded 



Runge-Kutta methods. Both methods in the pair are based 
on piecewise-polynomial approximants, which are used for 
error estimation and step size selection, to handle the neces- 
sary interpolations for delayed solution values, and to track 
derivative discontinuities that are propagated by the system 
delays, while the sixth-order method is used to perform the 
actual integration. 



As discussed in (IShampine and Thom pson. 



2006) 



dde_solver was designed to solve systems of DDEs with 
state-dependent delays in as "user-friendly" a fashion as 
possible, while at the same time retaining and extending 
the solution capabilities of dklag6. Our solver was also 
designed s o that usage approaches the convenience of the 



MATLAB d2007l) DDE solvers dde2 3 and ddesd. For 
example, storage management is handled automatically by 
the solver, thus relieving the user of the necessity to supply 
work arrays whose sizes are not known in advance. Several 
options are available for supplying necessary information 
about the problem and for dealing with its special character- 
istics. All options have carefully chosen defaults that can be 
changed by the user. These include the ability to supply vec- 
tors or functions to define the delays and the initial solution 
history, the ability to specify points corresponding to known 
derivative or solution discontinuities, tracking delay-induced 
derivative discontinuities, the ability to cope with small 
delays, the ability to handle state-dependent events {e.g., 
times at which it is desirable to make qualitative changes or 
parameter changes in the underlying system of DDEs), and 
the ability to solve so-called neutral DDEs, which contain 
delayed derivatives. The solver builds and returns a Fortran 
90 solution structure that may be used for various tasks, e.g., 
for plotting purposes. An interpolation module uses this 



1997). Both dde_solver and dklag6 implement a structure, for example, to perform additional interpolations 
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requested by the user. 

The numerical studies in this paper led us to incorporate 
several new options in dde_solver. By default, the solver 
retains the entire solution history. The numerical simulations 
in this paper, though, require very long integration intervals, 
leading to prohibitive storage requirements. We added, there- 
fore, an option to have the solver trim points from the solu- 
tion history queue that precede the largest delay. A related 
option was added to allow the user to provide a module to 
process solution information before it is discarded, thus re- 
taining the ability for user interpolation. 

The dde_solver with these added options is available 
at http://www.radford.edu/~thompson/ffddes/ In addition to 
the solver, a variety of example programs can also be found 
there; they may be used as convenient templates for other 
problems. This solver is a Fortran 90/95 compliant self- 
contained module with no restrictions on its use. In par- 
ticular, it is not compiler dependent and has been used suc- 
cessfully with most of the available F90/F95 compilers, in- 
cluding, for example, g95, Lahey LF90, Lahey-Fujitsu LF95, 
Salford FTN95, SUNf95, and Compaq. 



Matlab problem solving environment; dde2 3 is based on 
a (2,3) pair and is applicable to DDEs with constant de- 
lays. It is worth noting that the dde2 3 user interface led 
to many of the design decisions used in dde_solver. The 
solver ddesd incorporates novel methods based on control 
of the solution residual and is intended for systems with state- 
dependent delays. In addition to these solvers, two other 
popular and well-known tools incl ude the DDE-BIFTOOL 



package (Engelborg hs et al. 



age (lErmentrout 



2001) and the XPPAUT pack- 



2002b . Each of these packages contains a 



variety of tools that are useful for analyzing delayed dynam- 
ical systems. 
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